Quantum optimal control theory and dynamic 
coupling in the spin-boson model 
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A Markovian master equation describing the evolution of open quantum systems in the pres- 
ence of a time-dependent external field is derived within the Bloch-Redfield formalism. It leads to 
a system-bath interaction which depends on the control field. Optimal control theory is used to 
select control fields which allow accelerated or decelerated system relaxation, or suppression of relax- 
ation (dissipation) altogether, depending on the dynamics we impose on the quantum system. The 
control-dissipation correlation and the non-perturbative treatment of the control field are essential 
for reaching this goal. The optimal control problem is formulated within Pontryagin's minimum 
principle and the resulting optimal differential system is solved numerically. As an application, we 
study the dynamics of a spin-boson model in the strong coupling regime under the influence of an 
external control field. We show how trapping the system in unstable quantum states and transfer 
of population can be achieved by optimized control of the dissipative quantum system. We also 
used optimal control theory to find the driving field that generates the quantum Z-gate. In several 
cases studied, we find that the selected optimal field which reduces the purity loss significatively is a 
multi-component low-frequency field including higher harmonics, all of which lie below the phonon 
cutoff frequency. Finally, in the undriven case we present an analytic result for the Lamb shift at 
zero temperature. 

PACS numbers: 32.80.Qk, 03.65.Yz, 78.20.Bh, 78.67.-n 
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I. INTRODUCTION 

In the theory of quantum information and computa- 
tion, quantum coherence and entanglement are used as 
essential resources for efficient information processing 1 . 
However, the interaction of the quantum system with its 
environment eventually leads to a complete loss of the 
information initially stored in its quantum state. This 
phenomena, known as decohcrence, is regarded as a seri- 
ous obstacle to a successful implementation of quantum 
information processing. The question of how it is possi- 
ble to avoid the negative influence of this process is one of 
the most interesting issues in modern quantum mechan- 
ics. It not only concerns the area of quantum information 
and computation but many other fields of physics as well. 
The challenge is to preserve quantum coherence during 
a sufficiently long time needed for both storage and ma- 
nipulation of the quantum states in systems which are 
unavoidably exposed to the influence of their surround- 
ing environment. 

Over the last few years, a number of interesting 
schemes have been proposed to eliminate the undesirable 
effects of decoherence in open quantum systems, includ- 
ing decoherence free subspaces^, quantum error cor- 
rection codesi^, quantum feedback^ and mechanisms 
based on the unitary "bang-bang" pulses and their gener- 
alization, quantum dynamical decoupling 8 9 - i^iiiiiSii&ii. 
The key ingredient of dynamical decoupling is the con- 
tinuous disturbance of the system, which suppresses the 
system-environment interaction. It has been shown that, 
in the bang-bang control schemes, the decoherence of 
the system is effectively suppressed if the pulse rate is 



much higher than the decoherence rate due to the system- 
environment interaction. As already pointed out by Vi- 
ola and Lloyd 8 ^, the situation is similar to the so-called 
quantum Zeno effect which takes place in a system sub- 
ject to frequent measurements projecting it onto its ini- 
tial state: if the time interval between two projections 
is small enough the evolution of the system is nearly 
"frozen" . In a similar manner to the quantum Zeno ef- 
fect, a fast rate control freezes decoherence. Analysis and 
comparison of the effects of these different physical pro- 
cedures (bang-bang dynamic decoupling, coherence pro- 
tection by the quantum Zeno effect and continuous cou- 
pling) have been investigated in Refii&. Advances in de- 
coherence control using dynamical decoupling startegies 
is addressed in Refit 

The staring point of the decoupling techniques is the 
observation that even though one does not have access 
to the large number of uncontrollable degrees of freedom 
of the environment, it is still possible to interfere with its 
dynamics by inducing motions into the system which are 
at least as fast as the environment dynamics 11 . More- 
over, if one can establish an additional coupling to the 
system by means of an external control, there can be 
quantum interference between the two interactions. The 
degree and nature of quantum interference - constructive 
or destructive - can be controlled by adjustment of the 
control fields. 

In a simple-minded model of a dissipative quantum 
system, where the interference between the system-bath 
and system-control interaction is ignored or is irrele- 
vant only limited control can be achieved^. The sit- 
uation changes dramatically when interference between 
the system-environment and system-control interaction 
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can be used to control the effective system-environment 
coupling 18 i 20 i 21 i 22 i 23 i 24 i 25 i 26 i 27 i 28 i 29 . This effect of coher- 
ent control of "dissipation" is demonstrated here for the 
example of the driven spin-boson model in which a quan- 
tum two-level system (qubit) is modelled by a spin, the 
environmental heat bath by quantum oscillators, and the 
spin subjected to an external control field is coupled to 
each bath oscillator independently. Decoherence control 
of this model is formulated using optimal control which 
is mathematically a problem of functional optimization 
under dynamical constraints^Si 3 ^* 3 ^. Recently, we stud- 
ied the same model with a control field restricted to a 
monochromatic wave plane^ 9 .. The task was to find a set 
of three parameters namely the amplitude, the frequency 
and the phase using optimal control theory. Results were 
presented for control of the relative population of the spin 
system, i.e, the z-component of the Bloch vector. In the 
present paper, we show how this work can be extended 
to a control of all components of the Bloch vector simul- 
taneously, as well as to general control field shapes. 



II. QUANTUM MASTER EQUATION FOR 
DRIVEN OPEN SYSTEMS 

Consider a physical system S embedded in a dissipa- 
tive environment (B also referred to as the heat bath) 
and interacting with a time-dependent classical external 
field, i.e., the control field. The Hilbert space of the total 
system 7itot = 7~Ls ® 7~Lb is expressed as the tensor prod- 
uct of the system Hilbert space Tts and the environment 
Hilbert space TCb ■ The total Hamiltonian has the general 
form 

fltot = H s (t) + H B + H int , (2.1) 

where Hs(t) is the Hamiltonian of the system, Hb of 
the bath and H mt of their interaction that is responsible 
for decoherence. The operators Hs(t) and Hb act on 
Hs and Hb, respectively. The operator Hs(t) contains 
a time-dependent external field to control the quantum 
evolution of the system. 



The spin-boson model is a widely used model sys- 
tem. It can be mapped to a number of physical situ- 
ations 3 ^. In the theory of open quantum systems, the 
spin-boson model is actually one of the most popular 
models and has gained recent practical importance in 
the field of quantum computation 1 . A special variant of 
it, in which the inter-level coupling is absent, is known 
as the independent-boson model^i. These models have 
been used to study the role of the electron-phonon in- 
teraction in point defects and quantum dots, interact- 
ing many-body systems, magnetic molecules, bath as- 
sisted cooling of spins and a two level Josephson- Junction 
^^Si^LSi 3 ^ 4 ^. Its basic properties have been reviewed in 
the literature 41 -. 



The remainder of the paper is organized as follows: 
in the next section we present a pedagogical derivation 
of Born-Markov master equations for dissipative N-level 
systems in the presence of time-dependent external con- 
trol fields. The master equation is written as a set of 
Bloch-Redfield equations and a qualitative discussion of 
the influence of the control field on dissipation is given. 
This equation is the starting point for the derivation of 
the kinetic equation for the driven spin-boson model in 
the strong spin-boson coupling regime, as outlined in 
Sec. lIIII In Sec. lIVI within the Pontryagin minimum prin- 
ciple we formulate the optimum control problem in terms 
of the Bloch vector. The general cost functional and its 
gradient in case of arbitrary control field are given. We 
also present the numerical approach in form of the gradi- 
ent method. Our numerical results presented in Sec. IVl 
Summary and conclusions are given in Sec. IVII Some 
mathematical details are relegated to appendixes. 



A. Perturbation theory in the system-bath 
coupling 

We shall be interested in the time evolution of the re- 
duced density matrix of an open system, defined as 

Ps(t) = tr B { P tot(t)} ■ (2.2) 

where ptot is the total density matrix for both the system 
and the bath. Here and in the following trs denotes the 
partial trace over the open system's Hilbert space, while 
tiB denotes the partial trace over the degrees of freedom 
of the environment B. A number of different methods 
are used to derive an equation of motion for the reduced 
density matrb*^ 2 -. However, in any practical applications, 
two approximations are commonly invoked. The first is 
the initial factorization ansatz. Basically, one assumes 
that the interaction Hi nt is switched on at time t = 0. 
Prior to this, the system S and the environment B are 
uncorrelated so that the initial density matrix of the total 
system factorizes as the product of the system and bath 
contributions, that is, 

Ptot(0) = ps(0)®pb (2.3) 

with p B = tr s {ptot(0)} and p s {0) = tr B {ptot(0)}. The 
second approximation is the weak system-bath interac- 
tion limit in which the second-order perturbation theory 
is applicable. In the following, we shall present a deriva- 
tion of the master equation for ps- H lrA is assumed to 
be weak and will be treated perturbatively up to second 
order for its effects on ps under coherent external driving. 

Let us start with the equation of motion for the total 
density matrix, 

Ptot(t) = ~[H tot (t),p tot (t)] . (2.4) 
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The formal solution of the above equation can be ob- 
tained in the interaction picture 2 ^ with respect to the 
Hamiltonian H(t) = H s (t) + H B , 



pUt) = UHt,Q) Ptot (t)U(t,0) 



where the unitary operator reads 



U{t,t') = T\ exp 



cLtH{t) 



(2.5) 



(2.6) 



The time ordering T of the exponential in Eq. 1)2.6)1 is 
defined as the Taylor series with each term being time 
ordered. Using the fact that the operators Hsif) and 
H B commute 



[H s (t),H B }=0 
U(t, t') can be decomposed in the form 

U(t,t') = U s (t,t')®U B (t-t') 

with 



U s (t,t')=T exp 



cItH s (t) 



(2.7) 
(2.8) 

(2.9) 



being the propagator of the coherent system dynamics 
satisfying the Schrodinger equation 



ih^Us{t,t') = Hs{t)Us{t,t') 



subject to the initial condition 

U s {t',t')=X 

and 

i 



U B (t-t') =exp 



H B (t-t') 



(2.10) 
(2.11) 
(2.12) 



is the propagator describing the free evolution of the en- 
vironment. The equation of motion of p{ ot is then 



pL(t) 



[Hi(t), d*®] 



with H](t) defined as 

H I {t) = UHt,0)H int U{t,0). 



(2.13) 



(2.14) 



Hj(t) is the interaction Hamiltonian written in the in- 
teraction picture. In integral form, Eq. (|2.13|) can be 
rewritten as 

p( ot (t) = Ptot(0) - % -J dt'[ ff,(f J.z&tCtO] (2-15) 

Inserting Eq. (|2.15(l into (|2.13(l and taking the trace over 
the bath we find the equation of motion for the reduced 
density matrix of the physical system 

Ps® = dt'tx B {[Hj{t), [Hjit'lpUt')]]} ■ 

(2.16) 



We have assumed that (using Eq. I|2.3|l ) 

tT B {[Hi(t),pt ot {0)]} = tr fl {[Hi(t),ps{0)®pB]} = 

(2-17) 

This assumption has to be justified when considering a 
specific case. Equation 1)2.16)1 still contains the density 
matrix of the total system p( ot (t) on its right-hand side. 
In order to eliminate p{ ot (t) from the equation of motion, 
we assume that the interaction between the system and 
the environment is weak and perform the Born approxi- 
mation^ 2 -* 3 . 



pU*') = ps(0 



Pb 



(2.18) 



In practice, one usually assumes a thermal equilibrium 
for the environment, 



-I3H E 



PB 



tr e 



(2.19) 



where f3 — l/k B T with T the environment equilibrium 
temperature. Inserting the tensor product l|2.18|l in the 
exact equation of motion (|2.16l) . we obtain a closed 
integro-differential equation for the reduced density ma- 
trix p s (t) 

= ~jp J dt'tr B {[H I (t),[H I {t l ) ) p 1 s {t l )®p B \]}. 

(2.20) 

In order to further simplify the above equation we per- 
form the Markov approximation 23 ' 42,43 in which the in- 
tegrant Ps(t') is replaced by Pgit) 

Ps(t) = -^ J dl/ti B {[H J (t),[H I ^),p l s (t)»pB]]} ■ 

(2-21) 

The Markov approximation is therefore justified if the 
time scale tr over which the state of the system varies 
appreciably is large compared to the time scale t b over 
which the reservoir correlation functions decay (tr 3> 
t b ). One now goes back to the Schrodinger picture using 
Eqs. 1)2. 5[l . I)2.14[l . the decomposition l|2.8)l . and the cyclic 
invariance property of the trace 



ps(t) = --[H s (t),p s (t)} 



1 



dt'tr B { ^H int (t - t'), [H- m t{t, 0. Ps(t) 



Pb 



}■ 

(2.22) 



Here 



H int (t,t') = U s (t,t')H int Ul(t,t') 



(2.23) 



is the interaction Hamiltonian in the anti-Hciscnbcrg rep- 
resentation with respect to Hs(t) and 



fli* (t - = u f B {t- t')H int U B (t - 1') 



(2.24) 



the interaction Hamiltonian in the Hcisenberg represen- 
tation with respect to H B . The time-evolution operators 
U s (t,t') and U B (t - f') are defined by Eqs. and 
(|2.12|) . respectively. 
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B. Bloch-Redfield equations 

Let us suppose that Tts is an ^-dimensional Hilbert 
space of orthonormal basis \i),i = 1 ... N. Writing the 
quantum master equation i|2.22[l in the representation \i) 
and expanding the double commutator, we obtain after 
some algebra-2^ 

ps.ij(t) = ~T ^ (Hs,ik(t)Sij - 5ikHs,ij(t)) ps,ki(t) 

kl 

-^n ijU {t) p s ,kl(t) (2-25) 

kl 

where the first term represents the unitary part of the 
dynamics generated by the system Hamiltonian Hs(t) 
and the second term accounts for the dissipative effects 
due to the coupling of the system to the environment. 
The Redfield relaxation tensor Riju (t) is given by 

TZ ijkl (t) = S l3 ^2 r irrk(t)+^kY, r ir r] (t) 

r r 

-Itf ifc (*)-ry ifc (t) (2.26) 
with the time-dependent rates Tfj kl (t) 

(2.27a) 
(2.27b) 

Let us now write the Hamiltonian iJ; nt in the Schrodinger 
picture in the following general form 

H in t = ^A a ®B a (2.28) 

a 

where A a and B a are hermitian operators of the system 
and the environment, respectively. Inserting the expres- 
sion lf!T2"gjl into Eqs. (|2.27a[) and (|2.27bj) . using the fact 
that A a and B a commute leads to2M£ 

T hk( t ) = i [ dt'Y J (B a {t~t')B p (Q)) B A a , l] 

a,/3 

X ^ U S , im (t, t')Ap 

.ran. 

(2.29a) 

r y,*(*) = P / dt'Y J {B fJ {0)B a {t-t')) B 

.III 1 1 

m,n 

(2.29b) 



where 

(B a (r)B p (0)) B - tr B {B a ( T )B p {0)p B } (2.30) 

is the environment correlation function with 
(B a (0)B p (T)) B - (B a {-T)B p {0)) B , since the envi- 
ronment is supposed to be in thermal equilibrium and 
its evolution is then homogeneous in time. Note that 
the condition H2.17[) becomes 

(B a (r)) B =ti B {B a (T)p B } =0, (2.31) 

which states that the reservoir averages of B a (t) vanish. 

The Born-Markov master equation l|2.25|l with the 
time-dependent decay rates defined in Eqs. (j2.29a.fl and 
(I2.29b(l together with the Schrodinger equation (|2.10(l 
satisfied by the coherent time evolution operator Us(t, t') 
provide all the necessary ingredients to describe the dy- 
namics of a driven open quantum system. Note that the 
interaction of the system with the time-dependent con- 
trol Hamiltonian Hs(t) is treated non-perturbatively in 
the derivation of the above quantum master equation. 
The application of the present formulation to the driven 
spin-boson model in the weak coupling regime enables 
us to find an identical set of Bloch-Redfield equations 
obtained from projector-operator methods by Hartmann 
et alA 

In the absence of the control field, Us(t, t') = Us(t — t') 
is the free time evolution operator of the physical sys- 
tem. In this particular case, the time integration in Eqs. 
(|2.29af) and (|2.29bl) can be replaced by \ in the 

Markov approximation. After the extention of time in- 
tegration to infinity the decay rates T^ kl become time- 
independent and arc then given in terms of the Fourier 
transform of the product of the environment correlation 
functions and the system dissipative operators (the inter- 
action vertices), in agreement with the well-known results 
for undriven quantum open systems^. 



C. Control field effects 

The analysis of Eq. (|2.25(l shows that the presence of 
a time-dependent external field affects both the unitary 
evolution and the dissipative parts of the quantum mas- 
ter equation. The dissipative field influence can be inter- 
preted as a direct consequence of quantum interference 
between the system-bath interaction and a coupling of 
the system to the external field. In order to study the dis- 
sipative field influence, let us examine the transition and 
the dephasing rates in the secular approximation^^. In 
the zero field limit, the secular approximation assumes 
that the populations and the coherences are completely 
decoupled. It is valid if the typical time scale of the evo- 
lution of the system rs, defined by a typical value for 
the inverse of the frequency associated with the system's 
energy levels, is small compared to the relaxation time 
of the system, that is t$ If we suppose that the 
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secular approximation is still valid in the non-zero field 
case, the transition rates are defined as^SL^ 



Wait) = nasi® 
(*)- 



31,13 



3hi3 



(*), 



(2.32) 



with i ^ j. Then, from the property 



Ij.ik 



(*) 



it follows that the time dependent parameters 
Wy'(t) are real, W£(f) = Wy(f). On the other hand, the 
dcphasing rates are defined as22i2L^ 

r«(t) = = £ [r£ iri (i) + rT >r .(t)] 



f (*) 



,•(*)• 



(2.33) 



The hermiticity condition of the density matrix implies 
^ij (t) ~ (i) which means that Tij (t) are complex num- 
bers. In the zero field limit, the parameters Wij lead to 
a population relaxation into a thermal mixture of the 
system's energy eigenvalues. Therefore, the diagonal el- 
ements of the density matrix decay to the value given 
by the Boltzmann factors. The real part of the parame- 
ters Tij(t) describes the dephasing, namely the decay of 
the off-diagonal elements of the density matrix towards 
zero. Their imaginary part leads to a renormalization 
of the system Hamiltonian which is induced by the vac- 
uum fluctuation of the environment quantum operators 
(Lamb Shift). If a time-dependent external control field 
is applied, all these quantities become time-dependent 
(via the control field) , and an external control of the dis- 
sipation is then possible. In particular, the correlation 
between the control field and the dissipation leads to the 
destruction of the detailed balance 



lim 



Wij{t) , exp{-0Ei) 



"00 Wji(t) 



exp{-0Ej) 



(2.34) 



and allows for the control of the states populations. Here 
Ei are the energy eigenvalues of the undriven physical 
system. Eq. (|2.34|l shows that the steady state can be 
far from equilibrium in the presence of a control field. 

Taking the limit tb — > as a reasonable approxima- 
tion, gives 



(B a (T)B p (0)) B cx S a0 6(T) 



(2.35) 



A random interaction with a ^-correlation function is 
called white-noise, because the spectral distribution 
which is given by the Fourier transform of (|2.35() is then 
independent of the frequency 45 . Substituting Eq. (|2.35|) 
into Eqs. (|2.29a|) and il2.29b(l . any field dependence dis- 
appears because Us(t,t) —1. For classical problems the 
white-noise approximation holds when, in the high tem- 
perature limit, the environment resides within the classi- 
cal regime and quantum effects may be ignored. In such a 
situation, the control-dissipation correlation disappears. 



III. DRIVEN SPIN-BOSON MODEL 

A. The model 

The driven spin-boson model consists of a two-level 
system interacting with a thermal bath in the presence of 
a time-dependent external controlS&ii. The Hamiltonian 
for this model can be written as 

fltot = H(t) + flint = H S (t) + H B + flint , (3.1) 

where the dynamics of the system S is described by the 
Hamiltonian 

H s (t) =~ (Act, + e Qz a z ) - ^e z (t)a z . (3.2) 

Here a a with a = x,y, z are the Pauli spin matrices; 
HA is the tunnelling splitting, Hso z is an energy bias and 
He z (t) is its modulation by a time-dependent external 
control field. The Hamiltonian of the environment is as- 
sumed to be composed of harmonic oscillators with nat- 
ural frequencies LOi and masses rrn, 



H R = 



Pi 



N 

^ V 2)1!; 

l—l 



2 2 9 



(3.3) 



where (x±, xn ,Pi, ■■-,Pn) are the coordinates and their 
conjugate momenta. The interaction between the system 
S and the environment B is assumed to be bilinear, 



flir 



N 



(3.4) 



where Cj is the coupling constant between the spin coordi- 
nate and the ith environment oscillator with coordinate 
qi while go measures the distance between the left and 
right potential wells. The coupling constants enter the 
spectral density function J (to) of the environment defined 
by, 



2 ^ 



c, 



mi u>i 



■5(u) - Wi) . 



(3.5) 



B. Polaron transformation 

The evaluation of the time-dependent Bloch-Redfield 
tensor lZijki(t) for the Hamiltonian 13.1fl . requires knowl- 
edge of the propagator of the coherent system dynamics 
Us(t, t'). Obtaining an analytical expression for Us(t, t') 
is not trivial because the Hamiltonian of the physical 
system (|3.2p . is time-dependent and not diagonal. To 
get round this difficulty we perform the polaron trans- 
formation of the Hamiltonian l|3.1|l . This transformation 
is defined by the unitary operatos^i 



V 



(3.6) 
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with 



ft = ^fii 



fti = {qoCi/Hm % u>1 



Pi 



Applied to the original Hamiltonian l|3.1|) leads to 

HU = VH^V- 1 = H' s (t) +H' B + H! at . 
The expression 



(3.7) 



(3.8) 



(3.9) 



is the internal system part of the transformed Hamilto- 
nian 



4 = iyfi 

B 9^1 m, 



defines the Hamiltonian of the bath, and 



-in 



2 ) 



(3.10) 



(3.11) 



gives the modified interaction. Here a± = (a x ± ia y )/2. 
After applying the polaron transformation the coherent 
propagator 



U s (t,t')=T lexp 



drH' s (r) 



(3.12) 



simplifies to 



U s (t, t') = cos [{e 0z (t - + f(t, t')} /2}1 

+^ sin [{e 0z (t - + f(t, t')} /2] a z (3.13) 

where X is the unit matrix of order 2x2 and 

f(t,t')=J i dre z (T). (3.14) 

A constant phonon-induced energy shift in the system 
Hamiltonian has been omitted. As an application of the 
general formulation of the kinetic equation for driven 
open systems developed in Sec. (jHjl, we will consider 
the Hamiltonian Ij3.8[l and derive the explicit form of the 
corresponding master equation for small A. 

The physical situation described by this model can be 
envisioned as a weakly coupled (via A) semiconductor 
double quantum dot containing a single electron. Each 
quantum dot provides one electronic level (spin is ig- 
nored). A metallic gate is located under each quantum 
dot. The relative bias between the two electrodes rep- 
resents the control field. Alternatively, the Hamiltonian 
studied here can be interpreted as a spin 1 /2-particle in 
a weak static B-field (eo, A) which is directed slightly 
away from the z-direction in an external control field ap- 
plied in z-direction. 

For the transformed Hamiltonian, the interaction con- 
tribution is not the system-bath interaction, but rather 



accounts for a phonon-renormalized coupling between 
the two states "up" and "down". Here, the dissipative 
operators of the system and those for the environment 
are Si = hAa+/2, S 2 = S\ = hAa_/2 and B 1 = e~ in , 
B2 = B\ — e* n , respectively^. For the new interac- 
tion Hamiltonian, Eq. (|3.11() , the rates Yf- kl (t) defined 
by Eqs. I|2.29a|) and (j2.29b|) may be written in terms of 
the equilibrium correlation functions with respect to the 
bosonic bath spectral density J{lo) in Eq. 13. 511 . 

( e =un(t) e ±ifi(o)^ = e -« a (t) e «ji(t) ( (315a) 

( e ±«(*) e^(°)) B = e -0»(*) e -«3i(*) , (3.15b) 
( e ±iQ(o) e ±<n(t)) B = e -« a (t) e -«3i(t) f ( 315c ) 

(e ±ffi (°) e^ n «) B = e-^C'e^W, (3.15d) 

where the exponents are given by^ 



(3.16a) 



„2 i-oo t/ \ 

Q 2 (t) = *| / du-\r(l-cos{ujt))coth(huj/3/2). 



nli 



C. Ohmic spectrum of the bath 



(3.16b) 



Within the present work, we consider the case of 
Ohmic spectrum 



(3.17) 



Here 7/ is a phenomenological friction coefficient while a 
is the dimensionless coupling constant, and tu c is an ultra- 
violet frequency cutoff. For the Ohmic bath H3.17fl . the 
functions Qi(r) in Eq. I|3.16a|) and Q2(t) in Eq. (|3.16b|) 
can be determined explicitly^ 



[(r) = 2a arctan(w c r), 



(3.18a) 



h{r) = 2a In 



r 2 (i 



(i + izb)vT+^ 5 



r 1 



' tit) 



r 1 



"hp j 
(3.18b) 



where T is Euler's Gamma function. 

Let us determine the behaviour of the function Q2{t) 
at low temperature. Using the relation 



1 ' T \ T ( T \ 7rT ^ 

1 Hp) \ +l h(3) = ~f3K sinh (nr/ph) 



(3.19) 



one obtains from Eq. I|3.18b(l for huj c (3 ^> 1 (low temper- 
ature) 



»(r) = aln(l+^r 2 ) + 2o:ln 



(3H ( 7TT 



(3.20) 



The first term is independent of temperature and de- 
scribes how the fluctuations of the field vacuum affect 
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the coherence of the open system. This contribution of 
Q2(t) to dissipation depends on the cutoff frequency u> c . 
The second term is the thermal contribution of Q2 (t) to 
dissipation. Notice its dependence on the thermal corre- 
lation time tb = Ph/n. On the other hand, the function 
Qi(t) is independent of temperature. 



D. Master equation 

For the description of the dynamics of a two-level sys- 
tem, it is convenient to map the state density matrix 



onto the Bloch vector p(t) 
defined by 



(Px(t),p y (t),p z (t)) 



Poi(t) + pw(t) 
p(t) = Tr(<rp(t)) = I t(poi(t) - P10W) | , (3.21) 
Poo(t) - pu{t) 

where cr = (a x , a y ,a z ) is the vector composed of the 
three Pauli matrices. Within this notation, the states of 
a two-level system are parametrized by a 3-component 
vectors in the Bloch-sphere B := {p € R 3 ; ||p|| < l}. 

By combining the Redfield equation ((2.25|) with 
Eq. 1)3.21(1 , we obtain for the Bloch vector the inhomoge- 
neous linear equation of motion, 

p(t) = (£ + e(t)) x p(t) + x (p(i) x n) 

-7(*)p(*)+7o(*) (3-22) 



with £ = (0,0,-£ Oz ) T , e(t) 
(0,1, 0) T , 



(0,0,-£ 2 (t)) T , n 




7oW 




(3.23) 



and the relaxation matrix 





7 (t) 
7(t) 



(3.24) 



where 
£(*) = 

7W : 

and 

70 (t) = 



-Im[W 1 2 1 i2(*)+ ^12,21 (*)] 

-im [r+ )21 (t) + ra ljl2 (t) - r+ jl2 (*) - rr 2 , 12 (t)] 

(3.25a) 

Re[Wii,n(t)-Wii l a 3 (t)] 

2Re[r+ )21 (t)+r+ )12 (t)] (3.25b) 



-Re[^u,n(t) + ^n,22(t)] 

-2Re[r+ i21 (t)-r+ il2 (t)] 



(3.25c) 



are linear combinations of the Redfield tensor elements. 
Eqs. (|3.15a|) - (|3.15d|) . together with Eqs. I|2.29a|) and 



(|2.29b|) for the rates Tfj kl (t) , and the analytical expres- 
sion of the coherent propagator ((3.13(1 lead to 



g-QaM sin [ £02T + ffc t _ T )] 0OB[Qi(t)] 

(3.26) 

cos[£ 0z t + f(t,t- t)] cos[Qi(r)] 



7(t) = A 2 f dTe- Q2(T ^cos[ 
Jo 



To W 



dre- Q2(T) sin[£ 0z r + /(t,i 



(3.27) 
rJJsinfQ^r)] 
(3.28) 



where the function /(t, i') is given in Eq. 1(3.14(1 while the 
functions Qi{t) and Qi{t) are defined in Eqs. (|3. l(ja(> and 
(|3.16b|) . respectively. 

The control field dependence of the rates defined in 
Eqs. 13.26(1 . ((3.27(1 and the inhomogeneous term Eq. 
(pH£|) enters via the function /(*,*') in Eq. l(O^I . The 
quantity £(t) describes renormalization effects on the sys- 
tem Hamiltonian since it depends on the imaginary part 
of the Redfield tensor elements. It serves as an effec- 
tive local-control field correction acting on the system. 
The relaxation and dephasing processes are determined 
by the rate f{t). Note that the values for £(t), ~y(t) and 
7 (t) at the current time t depend on the control field 
e(t") for t" € [0,t]. 

The explicit equations of motion for the components 
of the Bloch vector reads 



Px(t) = (£oz + £z(t))p y (t) 



(3.29a) 



Py (t) = ~(e 0z + e z (t)+Z(t))pAt)-l(t)Py(t), 

(3.29b) 

p z (t) = -7(*)P*(*)+7b(*)- (3.29c) 



We would like to emphasize here that the decoupling 
of the populations from the coherences follows from the 
polaron transformation and the perturbative expansion 
up to second order in the tunnelling coupling A, and thus 
no secular approximation is required. 



E. Undriven case 

In order to illustrate the effects of the bath, namely the 
relaxation of the system and its energy renormalisation, 
we can analyse the master equation in the absence of the 
control field. The analytical expressions for the rates are 
worked out in the Appendixes. At zero temperature and 
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for a > 1/2, the decay rate follows as 



7(£o 2 > Q)| T=0 
j(eoz < 0)| T=0 



ttA 2 



2r(2a) \u>, 



£ 



2a— 1 -(eoz/uc) 



Qz 



which states that the process of absorption of phonons 
and its inverse, the process of emission of phonons, oc- 
cur with equal probability in thermal equilibrium and 
arises from the following quantum property of the ther- 



(3.30a) ma l equilibrium correlation function (e J °( l ' e iS1 (°))s 



ttA 2 



1 



2r(2a) \uj, 



2d 



J(l(t-ih/3) p iO(0)\ 



(3.30b) F. Limits of validity of the polaron transformation 



which agrees in leading order in eQ Z /co c with the result of 
Ref4i. Here T is the Euler's Gamma function. A similar 
expression holds for the inhomogeneous term 



lo(£oz > 0)| 



vrA 2 



T=0 



2r(2a) 



2d 



2a- 1 -(e ,/u 
fc 0z e 



(3.31a) 



lo(£oz < 0)| 



7rA 2 



T=0 



2T(2a) 



(3.31b) 



Next, we consider the effect of the bath on the energy 
splitting. By using the expression of £( £ 0z) from l|C7[l. 
we obtain 



^Oz x 



£dz x 



1 



1 - 



tteoz) 



£0z 



T=0 



A 2 sinh {eoz/^c) T(2a - 2) 



(e *M) r(2o) 



(3.32) 

The last equation which is valid when ^ <§; 1 and a > 1 
shows that £(eoz)/eoz is negative and constitute one of 
the principal result of this work. The effect of £(eoz) is 
the analogue of the Lamb shift, i.e. the renormalization 
of the level splitting in atoms due to the coupling to the 
electromagnetic radiations. 

In thermal equilibrium, the system density matrix can 
be represented in the localised eigenstates \R) and \L) of 
the position operator a z — (\R)(R\ — \L)(L\) as 



= /3Ke 0z cr z /2 



P : 



2cosh(/3fe z/2) 



(3.33) 



The equilibrium values of the Bloch-vector can be cal- 
culated from the density matrix, p st = Tr (crp), yielding 



Pst 



\PxstiPy B t>P z 






tanh(S/3e 0z /2) 



(3.34) 



Equation l|3.34|l corresponds to the stationary solution 
of the master equation (|3.22(l . From the rate equation 
(|3.29c|l . it follows that the decay rate 7(eoz) and the in- 
homogeneous term 7o(eoz) satisfies the detailed balance 
condition 



7o(gQz) 

l{eoz) 



tanh(fi,/3e z/2) , 



(3.35) 



In the undriven case, the model with Hamiltonian 
(|3.1() cannot be solved analytically, in general, and there 
are no reliable approximate methods which apply for 
a fixed (maybe weak) coupling to the bath and for all 
temperatures including the very low ones. For sym- 
metric tunnelling (eq z — 0), application of perturba- 
tive theory in the Ohmic bath coupling leads to a non- 
analytical temperature dependence for the renormalized 
tunnelling A r cx T a . At higher temperature there is a 
crossover from damped oscillations to over-damped mo- 
tion^yj-ii; with a relaxation rate that, in the weak cou- 
pling regime (a <C 1), decreases with increasing tem- 
perature, 7 cx T 2q_1 . The singular behaviour of the 
weak coupling series shows that perturbative theories 
break down at low temperature. On the other hand, 
the method of polaron transformation with the resulting 
Hamiltonian l|3.8|) is basically a perturbation theory in 
the tunnelling parameter A and is suitable for the strong 
coupling regime as we have shown in the last section. 
In fact, the combination of the polaron transformation 
with the second Born approximation is equivalent to a 
double path integral non-interacting blip approximation 
(NIBAjSl. 

Recently, Vorrath et al. 4S used the combination of the 
polaron transformation with the Markov-Born approxi- 
mation to derive a master equation for the generalisation 
of the undriven spin-boson model to spins greater than 
one-half in the strong coupling regime. They showed 
that this method is good enough if the parameters of 
the model, namely the temperature and the coupling, 
are limited to the case where the NIBA is valid. In the 
case of the driven spin boson model, the limits of the 
NIBA are discussed in Ref^&. 



IV. QUANTUM OPTIMAL CONTROL 
PROBLEM 



Let the time t be in the interval [0, tp] for fixed tF- 
The evolution of the state variable p(t) governed by the 
master equation l|3.22|l depends not only on the initial 
state p(0) = pi but also on the time-dependent control 
variable e(t). The idea now is to seek the optimal form 
of the the control field that allows for steering the Bloch 
vector from the given initial state p(0) = pi to a desired 
final state at a specified time tp- Typically, it is possible 
to define a cost functional incorporating the objective. 
Then, the goal of optimal control algorithms is to calcu- 
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late the control field which can induce a specific dynamics 
by minimizing this cost functional constraint by the state 
equation j.e, the master equation (j3.22(l subject to 
the initial condition p(0) = p/. 



A. Cost functional 

For the problems of interest here, the cost functional 
may be written as 



J = * [p(*f)] 



If 



C(p(t),e(t))dt. (4.1) 



The functionals $ [p(tp)] and £(p(i),e(i)) account for 
the specific objective at the fixed target time tp and at 
intermediate times t £ [0,ip], respectively. J in Eq. (|4.1|) 
is the sum of the so-called final time cost functional and 
running cost functional. 



B. Pontryagin's minimum principle 

Consider the quantum optimal control problem of min- 
imizing the cost functional l|4.1|l subject to the dynamical 
constraint (|3.22|) : 



min { J = $ [p(t F )} + f tF C (p(i), e(t)) dt] 



p(t) - (e + e(i)) x p(i) + £(f) x (p(<) x n) 
-7(*)p(*)+7o(*) 

p(0) = P/ , te[o,t F ] 



(4.2) 



An optimal solution of the problem 1)4. 2|l is characterized 
by first order optimality conditions in the form of the 
Pontryagin's minimum principle 19 ' 30,31 ' 32 . These condi- 
tions are formulated with the help of the Hamilton func- 
tion that has the following form in our problem: 

H{p{t),\(t),e{t))=C(p(t),e{t)) + 

A(t) . {(£ + e(t)) x p(t) + £(t) x (p(t) x n) 

-7(*)p(*) + 7o(*)}. (4-3) 

where A(t) is called the adjoint or the co-state variable, 
which is a Lagrange multiplier introduced to implement 
the constraint and thereby to render the variables p(t) 
and e(t) independent. Following a variation in analogy 
to Hamilton's variation principle of classical mechanics, 
the necessary first order optimality conditions result in a 



two-point boundary value problem: 



■ = dH( P (t),\(t),e(t)) 

V ' 6p(i) ' 



te [0,t F ] 



(4.4) 



p(0) = Pi, 



A(*f) 







3K(p(t),A(t),£(t)) 



9p(t P ) 



te [0,t F ] 



where the last condition is equivalent to the vanishing of 
the first variation of the cost functional, i.e., SJ = 0. 

The minimum principle requires the solution of com- 
plicated nonlinear algebraic equations, namely, the opti- 
mality condition dTi/de — 0, which can only be solved in 
an iterative manner. The present optimal control prob- 
lem H4.21 is not singular because d 2 H/de 2 ^ 0, since £(t), 
■y(t) and f (t) depend nonlinear ly on e{t) regardless of 
the choice of the running cost C Applying then the im- 
plicit function theorem one concludes that the optimality 
condition may have one solution, i.e, e(t) = (p(t),\(t)) 
or more solutions which are only " locally unique" . Here 
we apply the gradient method as an iterative scheme for 
solving H4.2|) . 

Let us now explicitly compute the gradient of the cost 
functional with respect to the control field. For this we 
first write the equation of motion for the adjoint state A. 
Eqs. {OJ and 1031 lead to 



A(t) = - 



dC (p(i),e(t)) 



+ (e + e(t)) x A(f) 



dp(t) 

-K(()xA(())xn + 7 T (()A((). (4.5) 
or in components form 
dC 



Ax(t) : 

Xy(t) = ^ 



dp x {t) 

dC 

d Py (t) 

dL 

dp,(t) 



{e 0z +e z (t)+i(t))X y {t), (4.6a) 
{e Qz + e z {t)) \ x {t) + 1 {t)\ y {t) , 



•7(*)A,(<). 



(4.6b) 
(4.6c) 



Within the Pontryagin's Minimum Principle, the vari- 
ation of the cost functional (|4.1|l reads 



SJ 



SJ 



Sedt 



Se jo 
dC (p(t),e(t)) 



de(t) 



&H(p(t),X(t),e(t)) 
de(t) 

■f (A(t) x p(t)) . Se(t)dt 



Se(t)dt 



d 



de(t) 

d 
de{t) 



[A(i) ■ m) x ( P (t) x n)}] . Se(t)dt 

[A(t).{-7(t)p(t)+7o(*)}]-*e(*)*- 

(4.7) 
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TABLE I: Best-fit parameters for a model defined by 
Eq. 15.51 . A is the unit of the amplitudes A it A/2ir of the 
frequencies Vi and rd of the phases fa. \ 2 — 36.43. 



Parameters 


Fit 


Error 




V 


0.354089322 


2.431734350x10" 


5 


A x 


2.38068007 


6.047635201x10" 


3 


<t>i 


-1.83291567 


2.978756342x10" 


3 


A 2 


2.79811789 


6.033930931x10" 


3 


fa 


0.960902574 


3.718809069x10" 


3 


A 3 


-0.909244535 


6.041502784x10" 


3 


fa 


0.65358378 


8.091472983x10" 


3 


A 4 


0.448104715 


6.048448725x10" 


3 


fa 


-5.91914473 


1.480647739x10" 


2 


A 5 


0.144629845 


6.054115543x10" 


3 


<t>5 


-9.41113677 


4.260009231x10" 


2 


A 6 


-0.0758863206 


6.060940693x10" 


3 


fa 


2.94649922 


8.029036597x10" 


2 


A 7 


-0.0170884169 


6.059845970x10" 


3 


fa 


-0.882323605 


3.542289951x10" 


1 



The last equation enables us to compute the gradient of 
the cost functional with respect to the control field. 

We may summarize the whole procedure for computing 
the gradient of the cost functional as follows: i) for a 
given control field e(t), we first solve the state equation 
(|3.22() forward in time, ii) the solution obtained for p 
is then used for the backward integration of the adjoint 
equation (|4.5|) . iii) with p and A obtained we compute 
the gradient. 



C. Discretization 

By an appropriate discretization scheme, the above 
infinite dimensional constraint optimal control problem 
can be transformed into a finite dimensional optimization 
approximation^^. For this purpose, we subdivide the 
time interval [t\ — 0,tp] using 

t j+1 =t j + At, j = 1,...,M -1 with At = tp/M 

(4.8) 

The values of the state, the adjoint and the control are 
evaluated at the mesh points tj 

(p, A, e) = (pi, . . . ,p M , Ai, . . . , A A /,£i, • • -,£m) T G K 9M 

(4.9) 

where the following notation p(tj) := Pj, A(tj) := Xj 
and e(tj) := £j is used. 

Adopting the Eulcr scheme for solving the state equa- 
tion 1|3.22[) and the adjoint equation (|4.5I) and by apply- 
ing the Riemann-rule integration to the cost functional 
(|4.1|) and to its variation l|4.7|l . we obtain the main tool 
for solving the time-discrete formulation of the quantum 
optimal control problem defined in Eq. I|4.2|l 



• state equation 

Pj+i = Pj + At [(e + £j) x Pj + ^ x ( Pj x n)] 



-At 7jPi+7oj 
for j = 1, . . . , M — 1 with pi = pj 
• adjoint equation 

dC {pj,£j) 



(4.10) 



Aj-i = \j - At 



dpj 



+ (e + £j) x A^ 



-At[(£ j x\ j )xn + >yj\ j ] 



(4.11) 



for j = M,...,2with X M = ^1. 



• cost functional 



A I 



J = ^(p M )+AtJ2C(p j ,£ j ) 
3=1 



(4.12) 



• variation of the cost functional 
M 



AJ = AtJ2 



3=1 
A I 



dei 



+At E ^:[^-{^ x (p. xn )}]- A ^ 



j,fc=i 



M d 

+At E TJi- A , - { ->.i^ • 



Ae k 



j,k=l 



• gradient of the cost functional 



(4.13) 



dJ_ 



At 



dC {pi,£i 
d£i 



\i x p, 



M d 

+ Ai E [ X 3 " {$3 x (Pj x n ) - W + Toj } 

(4.14) 



3 = 1 



fori= 1,...,M. 



The projection of (|4.14(l on the z-axis of the Bloch 
sphere leads to 



dJ 

de zi 



= At 



(Pxi ? Pyi 3 Pzi ) ^zi) 



de 

AI 

3 = 1 
M 

+At^ X zj 

3 = 1 



<9£ 7 - 97,- 



de Z i 



-Pzj 



de zi 



(4.15) 



for i = 1, . . . , M 

The matrices d^j/dei, d^/j/dei and d^oj/dei of size 
M x M are not diagonal because in the time-continuous 
problem £(£), "f(t) and 70 (t) at the current time t depend 
on the control field at t' < t via the function l|3.14|) . 
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D. Numerical method 




FIG. 1: Implementation of the Z-gate with pj = (1,0, 0) T 
and pa — ( — 1,0, 0) T . Depicted are the Bloch vector p = 
(px,p y ,Pz) T and the linear entropy S = | (l — ||p|| 2 ) as a 
function of time for undriven case and for driven by the opti- 
mal control field obtained by the conjugate gradient method, 
(a), (b) and (c) show, respectively, the results for p x , p y and 
p z while (d) show the results for S. The parameters used are 
a = 0.2, eoz = A, uj c = 20A and k B T = = HA. The final 
time is set as t_F = 20/A and the chosen time step is 10 _2 /A 
corresponding to M — 2 x 10 3 as the number of mesh points, 
i.e, the dimension of the optimal control problem. 



The set of equations needed to solve the optimal con- 
trol problem (|4.2|l are the discrete-time versions of the 
cost functional J{e z \ . . . ,e z m) denned in Eq. I|4.12fl . the 
equation of motion for the state and the adjoint variables 
given by Eqs. I|4.10J) and l|4.11J) . respectively, and the gra- 
dient of the cost functional VJ(e 2 i . . . , £ z m) in the form 
of Eq. PTTSy 

If we can compute the cost functional and its gradient 
at arbitrary points e z — (e z \ . . . ,s z m) t £ M M , the gen- 
eral form of the gradient algorithm for minimization is 
as follows^ 

1. Initialization: the initial guess e 1 € K M and the 
stopping tolerance tol > are given; set i=l. 

2. Stopping test: 

• integrate the state equation forward in time 
to find p. 

• integrate the adjoint equation backward in 
time to find A 

• compute the gradient VJ(e^); 
if |VJ(4)| < tol stop. 

3. Computing the direction: compute the descent di- 
rection cf e R M defined by VJ(4) ■ d? < 0. 

4. Line-search: find an appropriate stepsize fi z > 
satisfying J(e\ + fj? d l ) < J{e\) 



5. Loop: set e l z +1 
to 2. 



e\ + fi 1 d 1 ; increase i by 1 and go 
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FIG. 2: Implementation of the Z-gate with p/ = (1, 0, 0) T and 
Pd = ( — 1, 0, 0) T . The upper panel (a) shows the optimal con- 
trol field selected by the conjugate gradient method (CGM) 
vs. time while the lower panel (b) shows its power spectrum. 
A comparison with the model fit defined by Eq. (15.511 is also 
shown in the upper panel (a). Parameters as in Fig. Q 



In the work described here, the optimization of the 
cost functional is performed by using the subroutine FR- 
PRMN of the Numerical Recipes package^ which im- 
plements the conjugate gradient method as a variant of 
the above descent algorithm . We also used the subrou- 
tine DMNG of PORT library- 1 implementing the quasi- 
Newton method. These two iterative methods of opti- 
mization are very popular. Both of them require the 
gradient but differ in the calculation of the descent di- 
rection. 

The equations of motion for the state and the ad- 
joint variables are forward and backward initial value 
problems, respectively. We solved them using the Eu- 
ler scheme or a Runge-Kutta scheme which requires the 
values for the control field only at a grid point (see 
Sec. IIV C|) . Evaluation of the state and the adjoint 
variables involves an extensive computation of the time- 
dependent rates which are given by an integral over time 
of a rapidly oscillating functions (see Eqs. (|3.26l) . I|3.27p . 
(|3.28l) . H3.18al) and l)3.18bjl ). The numerical evaluation of 
the rate functions and their derivatives with respect to 
the control field involved in the computation of the gra- 
dient are performed using a Gauss quadrature suitable 
for an integration of rapidly oscillating functions. 
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FIG. 3: Implementation of the Z-gate with pi = (1,0, 0) T 
and pa = ( — 1,0, 0) T . The 3-dimensional plot of the the Bloch 
vector p = (px,p y ,Pz) T for undriven case and for driven by 
the optimal control field obtained by the conjugate gradient 
method is presented. Parameters as in Fig. Q 
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FIG. 4: Implementation of the Z-gate with pj = (1,0, 0) T 
and pd = ( — 1,0, 0) T . Depicted are the decay rate 7, the 
Lamb shift £ and the inhomogeneous term 70 as a function of 
time for undriven case and driven by the optimal control field 
selected by the conjugate gradient method, (a), (b) and (c) 
show, respectively, the results for 70, £ and 7. Parameters as 
in Fig. [I] 



V. NUMERICAL RESULTS 
A. Z gate 

As a first application of the quantum optimal control 
theory developed in Sec. I1VI we consider the action of 
the Z gate 




(5.1) 



which leaves |0) unchanged, and flips |1) to — 11). Its 
application to the initial state \ip) = (|0) + |1)) /y/2 leads 
to = (|0) — |1)) /v2- In term of the density matrix 
or the Bloch vector, we have for this particular state 





(5.2) 



The action of the dissipative Z-gate is phrased as an op- 
timization problem. At time tj = the two-level system 
(qubit) is prepared in the initial state pj = (1,0, 0) T . 
Our objective is to bring it into the desired state = 
(—1, 0, 0) T at time t = tp. In this case, we need to mini- 
mize the deviation of the state of the system at final time 
p{t F ) from the desired state pd- The cost functional cho- 
sen for this task is 



J 



Pd 



(5.3) 



corresponding to the running cost functional 
£(p(t),e(t)) = for all t G [0,t F ] and to the fi- 
nal cost functional $ [p(^_f)] = ^IIp^f) — Pd|| 2 in 
Eq. I|4.1|) . The cost functional defined in Eq. I|5.3j) 
requires A(t^) = p(^f) — Pd as the initial condition in 
Eq. (|4.5|) for the backward integration of the adjoint 
state variables A. 

Fig-Hshows the components of the Bloch vector versus 
time and the evolution of the linear entropy defined by 42 



S(*) = (l-||p(0|| 2 )/2. 



(5.4) 



The dashed lines give the result for the case of zero con- 
trol e(t) — 0. The solid lines give the results for the 
optimum field which was obtained by starting from a 
zero initial field and allowing 20 iterations. Fig. [21 shows 
the optimal field versus time, as well as its power spec- 
trum. It can be seen that the selected field performs 
several abrupt switch operations between initial state 
Pi = (1, 0, 0) T and target state pd — (—1, 0, 0) T to arrive 
at the target state at time tp. In principle the Z-gate 
operation is completed at approximately time t = 2.5A. 
However, here we are interested in preventing the de- 
crease of the Bloch vector over a prolonged period of 
time. The physical interpretation to the selected solution 
is the following: inspection of the kinetic equations for 
the Bloch vector Eas. l|3.29a|) . I|3.29bl and l|3.29aj) shows 
that a static field e ZO pt(^) = — £oz makes (p^,0,0) T , for 
— 1 < Px < +1 a stable ( "decoherence-free" ) subspace of 
the driven system. In this optimization run, the gradient 
selected a multiple switching version, whereby the system 
is, approximately, switched between the decoherence-free 
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FIG. 5: Implementation of the Z-gate with pj = (0,1, 0) T 
and pa — (0,-1, 0) T . Depicted are the Bloch vector p = 
(px,p y ,Pz) T and the linear entropy S = \ (l — ||p|| 2 ) as a 
function of time for undriven case and for driven by the opti- 
mal control field obtained by the conjugate gradient method, 
(a), (b) and (c) show, respectively, the results for p x , p y and 
p z while (d) show the results for 5*. Parameters as in Fig. 



states p(t) = (1,0, 0) T and p(t) = (-1,0, 0) T . Dissipa- 
tion essentially is initiated during the first switching op- 
eration when there is a small build-up in p y and p z com- 
ponent, as can be seen by inspection of Fig^d) showing 
the linear entropy of the system. The latter increases al- 
most linearly with time, however, at a greatly reduced 
rate when compared to the time-evolution of the un- 
driven system. The situation is complicated because p z 
has a thermal equilibrium state at around 0.46. As one 
sees in Fig. Q{c) > the optimal field succeeds repeatedly in 
driving the p z (t) back towards zero. The evolution of the 
3-d Bloch vector is shown in Fig. [3J While the control- 
free evolution rapidly spirals towards the thermal equi- 
librium state p s t = (0, 0, tanh(ft/3e z/2)) T the selected 
optimum control field is able to stabilize the Bloch vec- 
tor and eventually drive it very near to the target state 
Pd = (-1,0,0) T 

Fig|2Ia) displays the time evolution of the selected op- 
timal field. The repeated switching of the Bloch vec- 
tor is achieved by a nearly periodic field. The essence 
of the Z-gate operation is more or less contained in 
one period. The electric field oscillates about the value 
£opt(0) ~ — £oz to trap the system in state p = (±1, 0, 0). 
The switching is performed by a positive pulse which 
is optimized to rotate the Bloch vector into state p = 
(=Fl,0, 0) T . Then the field goes negative again to trap 
the system in this state. Performing more iterations will 
smoothen the oscillation about £o p t(0) ~ — £q z and re- 
duce the slope in the rise of the linear entropy. The 
analysis of Figs. ^ and [21 shows that the small oscillations 
of the control field about the value — Sq z between two 



FIG. 6: Implementation of the Z-gate with pj = (0,1, 0) T 
and pa = (0,-1, 0) T . The upper panel (a) shows the opti- 
mal control field selected by the conjugate gradient method 
(CGM) vs. time while the lower panel (b) shows its power 
spectrum. 



switching operations (two positive pulses) are reflected 
in the time evolution of the Bloch vector. 

The influence of the control-field on the dissipative 
part of the kinetic equations and the energy renormal- 
ization is displayed in Fig. 01 showing 70, 7, and £ ver- 
sus time for the driven and undriven case. The periodic 
structure of the optimal control field manifests itself in 
both of them. The renormalization term £ and 70 re- 
semble, essentially, a shifted and rescaled version of the 
control field itself. In this fashion they optimize support 
for the action of the electric field, in particular, when 
the latter rises to perform a switching operation. The 
minima of the relaxation rate 7, on the other hand, oc- 
cur when the control field becomes large. In this way, 
dissipation during the switching process is minimized. 

Fig Htb) displays the power spectrum of the selected 
optimal control field showing seven pronounced peaks at 
near equidistant frequencies. So, the selected optimal 
control field can be approximated by 



e z (t) = J2 An sin ( 27ra ^ + <t>n) 



(5.5) 



depending on 15 adjustable parameters which we deter- 
mine using a nonlinear least square method consisting of 
minimizing the % 2 merit function defined by 



1 M 

in 

1=1 



(5.6) 



where M — 2000 is the number of mesh points and 
e z (ti) is the optimal control field shown in Fig^a), 
solid line. The results are presented in Tab.||J The value 



14 



of the fit parameter v = (0.35409 ± 2.43173 x 10~ 5 ) A 
corresponding to the first peak of the power spectrum in 
Fig Gib). The remaining higher frequency peaks are lo- 
cated at about ni/, n = 2 . . .7. Tab. [I] shows that the am- 
plitudes A n satisfy \A 2 \ > \Ai\ > \A 3 \ > ... > \A 7 \ while 
the phases 4> n alternate in their sign. In Fig. Eta) we com- 
pare the optimal control field selected by the conjugate 
gradient method with the model defined by Eq. I|5.5|l . 

We also studied flipping from state pi = (0, 1,0) T to 
Pd = (0, — 1,0) T . The results are presented in Figs. 
IH1 The same picture emerges. The optimized field im- 
mediately drives the system into state p = (1,0, 0) T , 
performs switching between the decoherence-free states 
p = (1,0, 0) T and p = (— 1,0, 0) T , and finally transfers 
it into the target state p = (0,-1, 0) T . Actually, Fig El 
displays the time evolution of the selected optimal field 
and its power spectrum. It is seen in Fig. |HJa) that the 
optimal control field starts out positive value to trans- 
fer the system from pj = (0,1, 0) T to p = (1,0, 0) T 
and goes negative value (approximately — Sq z ) to trap 
the system in this state. The switching is performed by 
a positive pulse which is optimized to rotate the Bloch 
vector into state p = (— 1,0, 0) T . Then the field goes 
negative again to trap the system in this state. After, 
performing several abrupt switch operations between the 
free-decohrence states p = (1,0, 0) T , p = (— 1,0, 0) T , 
the control field value at the final time tj? is positive 
in order to transfer the system into the target state 
p = (0, —1, 0) T . Contrary to the first example, the con- 
figurations p = (0,p y ,0) T for — 1 < p y < +1 arc not 
stable under external driving by a negative static con- 
trol field s z (t) — — eoz- Thereby, the control optimum 
field value is positive at the beginning and also at end 
of the time evolution interval [0, tjr] allowing the trans- 
fer of the system from p = (0, 1, 0) T to p = (1,0, 0) T at 
the initial time and from p = (1, 0, 0) T to the target state 
p = (0,-1, 0) T at a final time as it is illustrated in Figs.[S] 
andEfa). FigHJb) shows that the power spectrum dis- 
plays seven pronounced peaks at equidistant frequencies 
similar to the first example. The fitting model defined 
by Eq. I|5.5|l can be used for this case too which it is not 
shown here for brevity. 

For the two examples of implementing a quantum Z- 
gate, The conjugate gradient method selects a "multi- 
component low-frequency" . This aspect of the optimal 
control field is remarkable. Firstly, the optimum field is 
a superposition of harmonics. This allows one to iden- 
tify rather small number of optimization parameters for 
a direct optimization scheme, such as a genetic code. 
Secondly, all essential frequency components lie below 
the the Ohmic cut-off frequency uj c = 20 A. Hence, we 
have shown that there are optimized solutions for decou- 
pling system from environment at lower frequence than 
required in the "bang-bang" approach. 

The presented solutions was obtained by starting from 
control field zero and the optimization algorithm ob- 
tained, within the specified cost functional, a solution 
which performs 7 switching operations. In principle, one 
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FIG. 7: Implementation of the Z-gate with pi = (0,1, 0) T 
and pd = (0, —1, 0) T . The 3-dimensional plot of the the Bloch 
vector p = {Px,p y ,Pz) T for undriven case and for driven by 
the optimal control field obtained by the conjugate gradient 
method is presented. Parameters as in Fig. Q 




FIG. 8: Implementation of the Z-gate with pi = (0,1, 0) T 
and pd = (0,-1, 0) T . Depicted are the decay rate 7, the 
Lamb shift £ and the inhomogeneous term 70 as a function of 
time for undriven case and driven by the optimal control field 
selected by the conjugate gradient method, (a), (b) and (c) 
show, respectively, the results for 70, £ and 7. Parameters as 
in Fig. □ 



switching operation would be sufficient. Due to the pos- 
sibility to dynamically create stable intermediate states 
one is in a similar position as with transferring an elec- 
tron in an isolated two level system. In the latter case, 
increasing the intensity of a resonant harmonic light field 
induces an increasing number of Rabi flip operations. 
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B. Trapping 

In the following two examples we study the control of 
the z-component of the Bloch vector, physically, corre- 
sponding to the spin direction or relative population of 
"up" and "down" states. First we consider trapping of 
the system in the excited state = (0, 0, 1). p x = p y = 
at the initial time ensures that the Bloch vector has van- 
ishing x- and y-component in the future, regardless of 
the control field applied. The problem becomes one- 
dimensional in the Bloch-vector space. The chosen cost 
functional is 

J=i^- I* dt\\p(t)~ Pd \\ 2 (5.7) 

Jo 

In this case, the running cost functional follows as 
C(p(t),e(t)) = \\p(t) — Pd\\ 2 and the final cost func- 
tional <& [p(iF)] is equal to zero. For the isolated two- 
level system there are several known ways of trapping 
a two-level system by an external control with a x — 
coupling. One can make the trapping state to the ground 
state of the system or one can apply a monochromatic 
high-frequency field with matched intensity to induce 
dynamic localization 52 . These strategies can be gener- 
alized and be applied to the dissipative two-level sys- 
tem2Sij£. B th strategies have in common that one tries 
to find a control field which makes the trapping-state 
to an element of the decoherence-free subspace of the 
driven system. Following the first strategy a static con- 
trol field can be found to make the state p = (0, 0,p z ) T 
for —1 < p z < +1 to the thermal equilibrium ground 
state of the driven system for given finite temperature. 
Alternatively, a high-frequency field can be used to dy- 
namically decouple the open quantum system from the 
bath. In the "bang-bang" method mentioned in the in- 
troduction, this is achieved with a control field whose fre- 
quency is (much) higher than the maximum frequency of 
the bath&2i2L2£. In the present model this is the phonon 
cut-off frequency lu c . Here we will show that a dynamic 
decoupling can be achieved by a field whose characteristic 
angular frequencies lie below uj c . 

In the present model, an oscillating control field leads 
to a rapidly oscillating integrand for 7(f) and 70 (t) lead- 
ing to small values for these two functions. Fig. EJd) 
shows the time evolution of p z . The dotted line shows the 
free evolution of the system into its thermal-equilibrium 
ground state within a time of about 20/ A. Starting 
from a guess for the control field in form of a Gaus- 
sian pulse, an optimized solution is obtained via the con- 
jugate gradient method which stabilizes the system in 
state p = (0, 0, 1) T rather well. Comparing, the ini- 
tial guess to the selected optimal control field one sees 
in Fig. OJa) that the oscillations of the Gaussian pulse 
get picked up and are amplified. In regions were the 
Gaussian factor suppressed the field the selected optimal 
field is less structured. Fig. shows the power spec- 
trum of original guess and the selected optimal control 
field. The main peak from the original guess gets am- 



plified and higher harmonics of the central frequency of 
the original guess are used to fine tune the control field. 
The selected field still shows clear features of the orig- 
inal guess. This is quite typical for solutions obtained 
within the conjugate gradient method when more than 
one solution exists. Figs. EJb) and (c) show that state 
trapping is indeed caused by dynamic decoupling in this 
case. 7(t) and jo(t) show high frequency oscillations of 
small amplitude about zero. 

To address the issue of convergence of the numerical 
procedure we show the cost functional l|5. 7fl versus the 
number of iterations in Fig. 1111 It can be seen that, 
starting from a mediocre guess, convergence is reached 
typically within 10 iterations. Moreover, convergence is 
strictly monotonic. Compared to direct approaches, such 
as a genetic code, this method requires significantly lower 
number of computations of the cost functional and sig- 
nificantly less computation time. Thus the investment 
of setting up the optimization scheme by introducing the 
co-state and the extra task of backward integration for 
the latter pays off in the end. Moreover, the present 
method makes feasible the selection of " arbitrary" opti- 
mized control fields, 1. e. an optimization of the control 
field at every mesh point in time. Due to the large num- 
ber of mesh points this would make a direct optimization 
approach computationally highly expensive. There are, 
however, non-linear programming approaches which may 
fair well for the present system—. 



C. Inversion of population 

As a final example we consider the task of driving the 
system from its thermal equilibrium state (0,0, p z ) into 
the pure "up" state (0,0, 1) and subsequent trapping in 
this state. The general cost functional given by Eq. 14. 1|) 
is adapted to the present task by setting 

*lp(tF)} = ^-\\p z (t f )- Pzd (t f W (5.8) 

and 

£(p(t),e(t)) = ^\\ Pz (t)-p zd (t)\\ 2 +\s(t)e 2 (t) . (5.9) 

s(i), wi and W2 with wi + W2 = 1 are real- valued weight 
factors to specify driving (wi — 1) and trapping (1V2 = 1). 
One can use the latter two to shift significance between 
driving into a target state and driving the system along 
a specified trajectory p z d(t). The function s(t) may be 
used to tailor the control pulse shape. In case of certain 
linear control problems the third term is necessary to 
make the problem regular 54 . — 1 < p z d{t) < +1 defines 
the "desired" trajectory of the system. For the present 
discussion we set p z d(t) = 1- 

Numerical results are shown in Fig. I|12|) . Let us first 
look at the undriven case for an initial state (0, 0, 1), dis- 
played by the dotted lines. It is seen in Fig. dJ)(d) that, 
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FIG. 9: Trapping the system in an unstable quantum state, 
i.e, p(t) = (1,0, 0) T for all t G [0,t F ]. Depicted are the control 
field e z , the decay rate 7, the inhomogeneous term 70 and the 
relative population p z as a function of time for three cases of 
undriven, driven by an harmonic field with a Gaussian shape 
(the guessed control field) and driven by the optimal control 
field selected by the conjugate gradient method, (a) shows the 
results for the control field s z , (b) and (c) show, respectively, 
the results for 7 and 70 while (d) show the results for p z . 
The parameters used are a — 0.2, Eoz = —A, uj c — 20 A 
and k B T = /3' 1 = HA. The final time is set as t F = 20/A 
and the chosen time step is 10~ 2 /A corresponding to M = 
2 x 10 3 mesh points, i.e, the dimension of the optimal control 
problem. 



on the time scale considered, there is rapid thcrmaliza- 
tion of the system into the equilibrium state at about 
(0, 0, —0.96). Except for oscillations at very short times, 
j(t) and 7o(t) are essentially constant in time. 

For the driven case we consider two situations. In the 
first we wish to prepare the system in the target state 
(0,0,1) at tp — 500 when the system initially is in the 
thermal equilibrium state (0,0,-0.964). We set W\ = 1 
and W2 = 0. Since the intrinsic time scale is faster than 
the target time there exist many solutions to achieve the 
task. Here we choose an initial guess in form of a har- 
monic field of low frequency (adiabatic solution) and op- 
timize this guess subsequently with the conjugate gradi- 
ent method using 300 mesh points for the control field. 
We use s(t) to suppress the control field at times around 
zero, as well as high intensities. Results for an optimal 
solution are shown by the solid lines in Fig. Ijl2(l . The se- 
lected optimal control fulfils the conditions imposed and 
leads to a gradual transition into the target state. In 
this particular case, the qubit-environment coupling has 
effectively been reduced over the undriven case. 

The second case considered is driving the system from 
its thermal equilibrium state into the target state (0, 0, 1) 
as fast as possible and subsequently trap it there. In the 
cost functional we set wi — and W2 = 1- Again a 




FIG. 10: Trapping the system in an unstable quantum state 
i.ep(t) = (1, 0, 0) T for all t G [0, tF]- Depicted is the compari- 
son of the power spectrum of the optimal control field selected 
by the conjugate gradient method with the power spectrum 
of the guessed control field. Parameters as in Fig. [5] 



low-frequency harmonic field is selected for the initial 
guess and s(t) is used to tailor the selected control field 
at times around zero, as well to limit its intensity. The 
results are shown by the dashed lines in Fig. (|12|l . The 
selected control field rises sharply from about zero to, es- 
sentially a plateau. j(t) and 7o(t) vary significantly only 
in the time during the transfer. Although high fields are 
suppresses around time zero, the selected optimal control 
field manages a more rapid transfer into the new equilib- 
rium case (dashed line in Fig. I|12fl (d)) than the undriven 
case (dotted line in Fig. I|12|) (d)). Hence, we show that 
we have been able to significantly increase the effective 
interaction strength. 



VI. SUMMARY AND CONCLUSIONS 

In this work we have presented dynamic control of 
open quantum systems as an optimization problem. The 
Bloch-Redfield approach was used to derive Marko- 
vian kinetic equations of a driven open quantum sys- 
tem whereby the external control was treated non- 
perturbatively. This approach leads to a Redfield tensor 
which accounts for a coupling between system and bath 
which contains a causal dependence upon the external 
control field. Indeed, the present approach is equivalent 
to the time-convolutionless projection operator method 
within second order in the system-environment cou- 
pling 42 . This control-field dependence of the effective 
system interaction allows steering of the open quantum 
system and its coupling to its environment beyond what 
is feasible within a semiclassical treatment of the environ- 
ment in which interference between the system-control- 
field interaction and the system-environment interaction 
is neglected^. 

This approach was applied to the spin-boson model in 
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FIG. 11: Trapping the system in an unstable quantum state 
i.e p(t) = (1,0, 0) T for all t £ [0,t F ]. Shown is the cost 
functional vs. the number of iterations. Parameters as in 
Fig. El 



the strong electron-boson coupling limit. Using the po- 
laron transformation, the kinetic equations for the Bloch 
vector were derived and analysed. They feature an effec- 
tive coupling in the spin system which is renormalized by 
the spin-phonon interaction and displays a causal depen- 
dence (non-local in time) on the control field. Analytic 
results for Lamb shift and decay time are presented for 
the zero temperature limit in the absence of the control. 
It is shown for several examples that both the stationary 
states of the driven open quantum system and the rates 
at which they are reached can be controlled to a large 
degree. 

Steering of the open quantum system is formulated and 
solved as an optimization problem via Pontryagin's min- 
imum principle which is based on the introduction of La- 
grangean multipliers in form of a co-state (adjoint state). 
The set of optimality conditions is solved iteratively us- 
ing a conjugate gradient method. Numerical examples 
show that it leads to a monotonic improvement in the 
cost functional. 

Several physical situations have been investigated nu- 
merically to demonstrate quantum-interference-based 
optimal control of open quantum system. The studies 
of a 7T rotation of the Bloch vector in the x-y plane and 
trapping along the z-axis have shown that optimal con- 
trol fields of moderate frequency (as compared to the 
phonon cut-off frequency) can be selected which signifi- 
cantly extend the lifetime of purity and, hence, improve 
the chance of successful completion of an error free quan- 
tum operation or the storage of a dissipative system in a 
fixed quantum state. The analysis of driving and subse- 
quent trapping into a quantum state, which for the un- 
driven system is highly unstable, at the example of pop- 
ulation transfer has shown that this task can be achieved 
by slowly varying fields for the present model. Moreover, 
the rate of transfer can be varied within limits set by the 
maximally obtainable effective coupling strength of the 
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FIG. 12: Control of relative population p z (t)'- (a) shows the 
selected control field, (b) shows 7(t), (c) shows 70 (t) and (d) 
shows the corresponding time evolution. The dotted lines 
are for e(t) —0 (undriven case), the solid line is for transfer 
from (0,0,-0.965) to (0,0,1), and the dashed lines are for 
transfer to and trapping in (0,0,1). The parameters used 
are a = 0.25, £ 0z = -2, A = 0.5, cu c = 10, t F = 500, and 
temperature /3 _1 =0.5. 



open quantum system. The latter is determined by the 
system, the environment, and the system-environment 
coupling. 

This analysis has also shown that the inverse problem 
of identification of optimal control fields in general has 
a large number of optimal solutions. Within the conju- 
gate gradient method the selected optimal solutions usu- 
ally show a remnance of the initial guess. The number 
of optimal solutions may be reduced by additional con- 
straints which may be used to select experimentally fea- 
sible solutions, such as fields with a gradual rise time, 
rather than abruptly turned on fields. In some cases, 
quite different fields can produce near equal results. For 
example, trapping in a quantum state may be obtained 
by applying a static control field which makes the trap- 
ping state to its (approximate) new ground state. In 
this case a decoupling of the system-environment inter- 
action is not necessary or even desirable. It is in fact the 
system-bath coupling which drives the system system 
into its new equilibrium state. As an important result 
this study has shown that state-specific optimal control 
can be achieved by time-dependent fields whose charac- 
teristic frequencies lie below the maximum characteris- 
tic bath frequency. Alternatively, high-frequency high- 
amplitude "bang-bang" control fields, reminiscent of the 
effect of dynamic localization, may induce dynamic de- 
coupling by making the effective coupling strength small. 

Optimization of a dissipative quantum gate poses a 
more complicated problem than the one addressed here 
since optimization should occur independent of the input 
stat o 56 ! 57 . Moreover, the output state (target state), in 
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general, depends on the input state. We find that an opti- 
mal control field critically depends on the input state. A 
bang-bang solution (which is probably difficult to imple- 
ment in experiment) can be envisioned whereby a high- 
frequency high intensity field is applied to suppress the 
effective system-environment coupling. However, such 
a field usually also has a direct coupling channel to the 
system which may cause problems when the control field 
is not perfectly suitable for the input state. Whether 
the present approach which is based on specific trajecto- 
ries can be extended to optimize quantum gates by some 
averaging procedure or whether an optimization should 
directly be aimed at the superoperator responsible for 
the time-evolution will require further investigation. 
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APPENDIX A: SCHWINGER AND FEYNMAN 
REPRESENTATIONS 



FIG. 13: The Lamb shift for the undriven case and zero tem- 
perature as a function of the Ohmic cutoff frequency u) a . The 
test of our analytical result, Eq. 1C7L by comparison of a 
direct numerical integration of Eq. (IC10 is shown. The pa- 
rameters used are the coupling a = 1.2, the tunnel amplitude 

560 MHz. 



635 MHz and the energy bias Sia 



The Schwinger and Feynman representations^ will 
play an important role in the determination of the de- 
cay rate and the Lamb shift (see below). 

The Schwinger representation involves the Euler 
Gamma function defined by 



T(u) = dt 
Jo 



Rc v > 0, 



(Al) 



Making the variable change Du — t in the definition of 
the Euler gamma function (|A1|I , leads to 



1 



T{v) 



duu v - l e- uD < 



Rei/>0 ReL>>0. 

(A2) 

The identity l|A2(l allows to write the denominators D of 
the propagator in form of an integral on the Schwinger 
parameter u. 

On the other hand the Feynman representation^ in- 
troduces a parameter x (Feynman parameter) to squeeze 
the denominator factors into a single polynomial form 
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1 f 1 

= / dx [Ax + B{l~x)\ 

■ B J Q 



(A3) 



APPENDIX B: DERIVATION OF THE DECAY 
AT ZERO TEMPERATURE 

Here we derive an analytical expression of the decay 
rates 



70o*) 



dr e~ Q2 (t) cos [Q 1 (t)] cos [e 0z r] 



7/(>o*) +lb{e 0z ) ■ 



(Bl) 



with 



A 2 

7/(^0*) = ^-Rc 



dr i 



\ (B2) 



and 



A 2 f°° 

7b(eo,) = — Re/ dr e -«="We <0l We fao « T (B3) 
2 Jo 

are, respectively, the forward and the backward decay 
rates. Note that 7b(eoz) — 7 / (~ £ 0z)- Substituting 
Eqs. (|3.18a|l and (|3~2l7|) into Eq. JB2J, we obtain 



7/(eo«) = . _ 

./o (1 — vjJ c t) 



dr. 



(B4) 



Now with the help of the Schwinger identity i|A2() . 
Eq. l|B4p can be written as 



7/Ooz) = 



A 2 



-Re / duu 2a - l e- u 



(B5) 



2r(2a) 
Using the fact that 

[ dTe- l{eo *- u = ir5(e 0z ~uuj c )-iPP ( 

Jo \e 0z -uuj, 

(B6) 

where the first term is the Dirac distribution and the sec- 
ond term PP denotes the Cauchy principal part of the 
integral J Q dr/ (sq z — ulo c ) and introducing the Heavi- 
side distribution 0(u) to extend the bounds of integration 
from — oo to +oo, Eq. (|B5p becomes 



lf( e 0z) 



nA 2 



2T(2a) w c 



+ 00 



duu 201 - 1 e~ u 6{u)5 [ u- — 



U c J 

(B7) 
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Evaluating the convolution product l|B7|l . one ends up 
with the following formula: 



7/(eoz) 



7rA 2 / 1 



2n 



The Heaviside distribution (or step function), insures 
that at zero temperature absorption of energy from the 
environment is not possible. The final result for the decay 
rate is then for a > |, 



7(eo 2 > 0) 



-A- / 1 • 2 



2T(2a) \u c 



ttA 2 / 1 • -" 



(B9a) 
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For the inhomogeneous term 



7o (£ 0z ) = A 2 / dr e~ Q ^ sin[Qi(r)] sin^r] 



(BIO) 

similar calculation of the decay rate leads for a > 1/2 to 
ttA 2 / I 



7o(eo z > 0) = 
7o(eoz < 0) = 



2T(2a) 
ttA 2 



2r(2a) 



(J^j e 2 £~ x e-C^M) (Blla) 
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At zero temperature, the detailed balance condition takes 
the following form 



7o(£qz) = I 1 if e 02 > 
7(eo 2 ) I -1 if e 0z < 



(B12) 
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FIG. 14: Implementation of the Z-gate with p/ = (1,0, 0) T 
and pi? = (— 1, 0, 0) T . The upper panel (b) shows the numeric 
comparison between Pontryagin's minimum principle (PMP) 
and the finite difference (FD) approximation for computing 
the gradient of the final cost functional in Eq. (15.311 evaluated 
at the control field displayed in the lower panel (b). The 
parameters used are a = 0.2, eo z = A, uj c = 20A, fcsT = 
/3~1 — hA. The final time is set as tF = 10/A and the chosen 
time step is 5 x 10 _2 A corresponding to M = 2 x 10 2 mesh 
points, i.e, the dimension of the optimal control problem. 



where u is the dimensionless Schwinger parameter. The 
last integral can not be computed using the residues the- 
orem since it is singular at u = and at u — when 
a < 1/2. 

Nevertheless the application of the Feynman identity 
[3j to Eq. iJOHl, leads to 



£{saz) = ~ Jt n £ °o / dxl[x,a 
T(2a) u 2 c 7 V 



(C3) 



with 



APPENDIX C: DERIVATION OF THE LAMB 
SHIFT AT ZERO TEMPERATURE 



Let us now compute the Lamb shift given by 

/>OC 

£(e 0z ) = A 2 / dre~ Q ^ cos[Qi(r)] sin[e 0z r] 
Jo 

(t) | e iQi(T) e ie «T _|_ e -iQi(T) e is zT\ 



A 2 

= — Im / (In 





I'OC 








Jo 



du- 



(C4) 



which after the change variable v = u + {2x — l)^j* is 
transformed to 



(CI) 



Substituting Eqs. I|3.18a|) and [j!T2"0"jl into Eq. (EB and 
using l|B6|) . we obtain 



1 1 x, a, ) = / (#u ( v — (2a; — 1) — 



Now, the approximation — <C 1 leads to 
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Using Eq. I|C6|I. we get from (|C3|) for a > 1 



&0z 



sinh (e 0z /uJc) T(2a - 2) 



(eoz/w c ) T(2a) 



(C7) 



Combining again the Schwinger identity (jA2|l with 
Eq. i|C6(l and after some algebra we obtain for i < a < 1 



2 w a 
x lim ( 1 



2a-2 
(2 Q -1) log(2)7-l) 



r(2 - 2a) 



2a 



1 

(C8) 



A such limit in the last equation does not exist. 

I summary, our prediction for the renormalization of 
the energy bias due to the Lamb shift at zero temper- 
ature; in leading order in <C 1) and in strong 
coupling regime (a > 1), is the following: 



t()2 



: £0z X 



1 



A 2 sinh (e Qz /uj c ) T(2a - 2) 



w c 2 (e 0z /uJc) T(2a) 



(C9) 



The agreement of the analytical expression for the 
Lamb shift, Eq. (|C7|) . with the numerical integration of 
Eq. ijClfl by Gauss quadrature is shown in Fig. fPjl) . 



APPENDIX D: NUMERICAL TEST OF THE 
GRADIENT 



In order to test the method of Pontryagin's minimum 
principle (PMP) given by Eq. (|4.15|) to compute the gra- 
dient of the cost functional, Eq. (|5.3(l , we compare it with 
the finite difference approximation (FD). Fig. 1141 shows 
the result of the gradient for a control field of the form 
e z (t) = A(t) cos(flt + <fi) with frequency fl, phase </>, and 
a Gaussian envelope A(t). We can see in Fig. 1141 that 
the error, ERROR = |V PMP J - V FD j| , is roughly zero 
except at the switching times of the control field where 
its amplitude is suddenly increased. The good agree- 
ment observed for this case occurs because the frequency 
fi = 10 _1 xAof the control field is low and causes a slow 
variation of the cost functional. In this case the finite 
difference approximation is numerically stable and gives 
good results compared to the adjoint-state method. In 
case of high control field, this good agreement is lost be- 
cause the cost functional varies very rapidly and renders 
the finite difference method numerically unstable. 
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